This R Markdown script assesses models for the PAL (probabilistic associative learning) task of the EMBA project.
# number of simulations
nsim = 250
# set the seed
set.seed(2468)
The following packages are used in this RMarkdown file:
## [1] "R version 4.5.1 (2025-06-13)"
## [1] "knitr version 1.50"
## [1] "ggplot2 version 3.5.2"
## [1] "brms version 2.22.0"
## [1] "bridgesampling version 1.1.2"
## [1] "tidyverse version 2.0.0"
## [1] "ggpubr version 0.6.1"
## [1] "ggrain version 0.0.4"
## [1] "bayesplot version 1.13.0"
## [1] "SBC version 0.3.0.9000"
## [1] "rstatix version 0.7.2"
## [1] "effectsize version 1.0.1"
## [1] "BayesFactor version 0.9.12.4.7"
## [1] "bayestestR version 0.17.0"
We perform prior predictive checks as proposed in Schad, Betancourt and Vasishth (2020) using the SBC package. To do so, we create simulated datasets where parameters are simulated from the priors. These parameters are used to create one fake dataset. Both the true underlying parameters and the simulated values are saved. Then, we create graphs showing the prior predictive distribution of the simulated values to check whether our priors fit our general expectations about the data. Next, we perform checks of computational faithfulness and model sensitivity as proposed by Schad, Betancourt and Vasishth (2020) and implemented in the SBC package. We create models for each of the simulated datasets. Last, we calculate performance metrics for each of these models, focusing on the population-level parameters.
These evaluations were performed for models with two groups; however, they should still apply to the same models with three groups.
First, we create the data. Then, all predictors are set to sum contrasts.
# load the trial order
df.trl = read_csv('data/PAL_scheme.csv', show_col_types = F) %>%
select(trl, emo, tone, difficulty, phase, expected)
# create the data
df.pal = data.frame()
groups = c("A", "B")
nos = 22 # number of subjects per group
for (g in groups) {
for (i in 1:nos) {
df.pal = rbind(df.pal,
df.trl %>% mutate(subID = sprintf("%s%02d", g, i),
diagnosis = g))
}
}
# add bogus rts
df.pal = df.pal %>% mutate(rt.cor = 1)
# drop the neutral condition for the analysis
df.pal = df.pal %>%
filter(expected != "neutral") %>% droplevels() %>%
mutate_if(is.character, as.factor) %>%
ungroup()
df.var = df.pal %>%
group_by(subID, diagnosis, phase, expected) %>%
summarise(
rt.var = sd(rt.cor, na.rm = T) + 1
) %>%
filter(expected != "neutral") %>% droplevels() %>%
mutate_if(is.character, as.factor) %>%
ungroup()
# set and print the contrasts
contrasts(df.pal$diagnosis) = contr.sum(2)
contrasts(df.pal$diagnosis)
## [,1]
## A 1
## B -1
contrasts(df.pal$expected) = contr.sum(2)
contrasts(df.pal$expected)
## [,1]
## expected 1
## unexpected -1
contrasts(df.pal$phase) = contr.sum(3)
contrasts(df.pal$phase)
## [,1] [,2]
## postvolatile 1 0
## prevolatile 0 1
## volatile -1 -1
contrasts(df.pal$difficulty) = contr.sum(3)
contrasts(df.pal$difficulty)
## [,1] [,2]
## difficult 1 0
## easy 0 1
## medium -1 -1
contrasts(df.var$diagnosis) = contr.sum(2)
contrasts(df.var$diagnosis)
## [,1]
## A 1
## B -1
contrasts(df.var$expected) = contr.sum(2)
contrasts(df.var$expected)
## [,1]
## expected 1
## unexpected -1
contrasts(df.var$phase) = contr.sum(3)
contrasts(df.var$phase)
## [,1] [,2]
## postvolatile 1 0
## prevolatile 0 1
## volatile -1 -1
In the preregistration, we noted the following population-level effects for the model of the reaction time variances: group, expectancy, phase and difficulty. However, the posterior fit for this model was suboptimal, therefore, we dropped the predictor difficulty.
# figure out slopes for subjects
kable(head(df.var %>% count(subID, expected)))
| subID | expected | n |
|---|---|---|
| A01 | expected | 3 |
| A01 | unexpected | 3 |
| A02 | expected | 3 |
| A02 | unexpected | 3 |
| A03 | expected | 3 |
| A03 | unexpected | 3 |
kable(head(df.var %>% count(subID, phase)))
| subID | phase | n |
|---|---|---|
| A01 | postvolatile | 2 |
| A01 | prevolatile | 2 |
| A01 | volatile | 2 |
| A02 | postvolatile | 2 |
| A02 | prevolatile | 2 |
| A02 | volatile | 2 |
kable(head(df.var %>% count(subID, expected, phase)))
| subID | expected | phase | n |
|---|---|---|---|
| A01 | expected | postvolatile | 1 |
| A01 | expected | prevolatile | 1 |
| A01 | expected | volatile | 1 |
| A01 | unexpected | postvolatile | 1 |
| A01 | unexpected | prevolatile | 1 |
| A01 | unexpected | volatile | 1 |
# code for filenames
code = "PAL-var"
# model formula
f.var = brms::bf(rt.var ~ diagnosis * expected * phase +
(expected + phase | subID)
)
# set weakly informative priors
priors = c(
prior(normal(4.50, 0.50), class = Intercept),
prior(normal(0, 0.50), class = sigma),
prior(normal(0.15, 0.15), class = sd),
prior(normal(0, 0.25), class = b)
)
# check if the SBC already exists
if (file.exists(file.path(cache_dir, sprintf("df_res_%s.rds", code)))) {
# load in the results of the SBC
df.results = readRDS(file.path(cache_dir, sprintf("df_res_%s.rds", code)))
df.backend = readRDS(file.path(cache_dir, sprintf("df_div_%s.rds", code)))
dat = readRDS(file.path(cache_dir, sprintf("dat_%s.rds", code)))
} else {
# stimulate some data
set.seed(2486)
gen = SBC_generator_brms(f.var, data = df.var, prior = priors,
thin = 50, warmup = 10000, refresh = 2000,
generate_lp = TRUE, family = lognormal, init = 0.1)
bck = SBC_backend_brms_from_generator(gen, chains = 4, thin = 1,
warmup = 1500, iter = 6000)
dat = generate_datasets(gen, nsim)
saveRDS(dat, file.path(cache_dir, sprintf("dat_%s.rds", code)))
# perform the SBC
res = compute_SBC(dat, bck, keep_fits = F,
cache_mode = "results",
cache_location = file.path(cache_dir, sprintf("res_%s", code)))
df.results = res$stats
df.backend = res$backend_diagnostics
saveRDS(df.results, file = file.path(cache_dir,
paste0("df_res_", code, ".rds")))
saveRDS(df.backend, file = file.path(cache_dir,
paste0("df_div_", code, ".rds")))
}
We start by investigating the rhats and the number of divergent samples. This shows that only 0 of 250 simulations had at least one parameter that had an rhat of at least 1.05 but 14 models had divergent samples (mean = 5.93). This suggests that this model performs well, but we will check for divergence issues in the final models.
Next, we can plot the simulated values to perform prior predictive checks.
# create a matrix out of generated data
dvfakemat = matrix(NA, nrow(dat[['generated']][[1]]), length(dat[['generated']]))
for (i in 1:length(dat[['generated']])) {
dvfakemat[,i] = dat[['generated']][[i]][[1]]
}
# compute one histogram per simulated data-set
dvfakematH = dvfakemat
dvfakematH[dvfakematH > 1000] = 1000
breaks = seq(0, max(dvfakematH, na.rm=T), length.out = 101)
binwidth = breaks[2] - breaks[1]
histmat = matrix(NA, ncol = dim(dvfakemat)[2], nrow = length(breaks)-1)
for (i in 1:dim(dvfakemat)[2]) {
histmat[,i] = hist(dvfakematH[,i], breaks = breaks, plot = F)$counts
}
# for each bin, compute quantiles across histograms
probs = seq(0.1, 0.9, 0.1)
quantmat = as.data.frame(matrix(NA, nrow=dim(histmat)[1], ncol = length(probs)))
names(quantmat) = paste0("p", probs)
for (i in 1:dim(histmat)[1]) {
quantmat[i,] = quantile(histmat[i,], p = probs)
}
quantmat$x = breaks[2:length(breaks)] - binwidth/2 # add bin mean
p1 = ggplot(data = quantmat, aes(x = x)) +
geom_ribbon(aes(ymax = p0.9, ymin = p0.1), fill = c_light) +
geom_ribbon(aes(ymax = p0.8, ymin = p0.2), fill = c_light_highlight) +
geom_ribbon(aes(ymax = p0.7, ymin = p0.3), fill = c_mid) +
geom_ribbon(aes(ymax = p0.6, ymin = p0.4), fill = c_mid_highlight) +
geom_line(aes(y = p0.5), colour = c_dark, linewidth = 1) +
labs(title = "Distribution of simulated values", y = "", x = "") +
theme_bw()
tmpM = apply(dvfakemat, 2, mean) # mean
tmpSD = apply(dvfakemat, 2, sd)
p2 = ggplot() +
stat_bin(aes(x = tmpM), fill = c_dark) +
labs(x = "Mean", title = "Means of simulated values") +
theme_bw()
p3 = ggplot() +
stat_bin(aes(x = tmpSD), fill = c_dark) +
labs(x = "SD", title = "SDs of simulated values") +
theme_bw()
p = ggarrange(p1,
ggarrange(p2, p3, ncol = 2, labels = c("B", "C")),
nrow = 2, labels = "A")
annotate_figure(p, top = text_grob("Prior predictive checks", face = "bold", size = 14))
Subfigure A shows the distribution of the simulated data with bluer bands being more likely than greyer bands. It shows a general distribution that fits our expectations, even though there are a few values that are unreastically large. The same applies to the distribution of the means and standard deviations in the simulated datasets in Subfigure B. We go ahead with these priors and check the results of the SBC. We only plot the results from the models that had no divergence issues.
# get simulation numbers with issues
mx_rnk = max(df.results$max_rank)
check = merge(df.results %>%
group_by(sim_id) %>%
summarise(rhat = max(rhat, na.rm = T),
max_rank = mean(max_rank)) %>%
filter(rhat >= 1.05 | max_rank != mx_rnk),
df.backend %>% filter(n_divergent > 0), all = T)
# plot SBC with functions from the SBC package focusing on population-level parameters
a = 0.5
df.results.b = df.results %>%
filter(substr(variable, 1, 2) == "b_") %>%
filter(!(sim_id %in% check$sim_id))
p1 = plot_ecdf_diff(df.results.b) + theme_bw() + theme(legend.position = "none") +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 3))
p2 = plot_rank_hist(df.results.b, bins = 20) + theme_bw() +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 3))
p = ggarrange(p1, p2, labels = "AUTO", ncol = 1, nrow = 2)
annotate_figure(p, top = text_grob("Computational faithfulness and model sensitivity",
face = "bold", size = 14))
Next, we check the ranks of the parameters. If the model is unbiased, these should be uniformly distributed (Schad, Betancourt and Vasishth, 2020). The sample empirical cumulative distribution function (ECDF) shows deviations from the theoretical distribution (95%) and the rank histogram also shows ranks outside the 95% expected range. However, we judge this to be acceptable as we cannot identify clear bias patterns as described in the information accompanying the SBC package: https://hyunjimoon.github.io/SBC/articles/rank_visualizations.html
p3 = plot_sim_estimated(df.results.b, alpha = .8) + theme_bw() +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 3))
p4 = plot_contraction(df.results.b,
prior_sd =
setNames(c(0.50,
rep(0.25,
length(unique(df.results.b$variable))-1)),
unique(df.results.b$variable))) +
theme_bw() +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 3))
p = ggarrange(p3, p4, labels = "AUTO", ncol = 1, nrow = 2)
annotate_figure(p, top = text_grob("Computational faithfulness and model sensitivity",
face = "bold", size = 14))
Then, we investigated the relationship between the simulated true parameters and the posterior estimates. Although there are individual values diverging from the expected pattern, most parameters were recovered successfully within an uncertainty interval of alpha = 0.05.
Last, we explore the z-score and the posterior contraction of our population-level predictors. The z-score “determines the distance of the posterior mean from the true simulating parameter”, while the posterior contraction “estimates how much prior uncertainty is reduced in the posterior estimation” (Schad, Betancourt and Vasisth, 2020). The contraction reveals very narrow posterior standard deviations with slightly lower contraction scores than we would want in the optimal case, however, we judge this as acceptable.
In the preregistration, we noted the following population-level effects for the model of the reaction time variances: group, expectancy, phase and difficulty; as well as the group-level predictores subject and trials.
# figure out slopes for subjects
kable(head(df.pal %>% count(subID, expected)))
| subID | expected | n |
|---|---|---|
| A01 | expected | 240 |
| A01 | unexpected | 48 |
| A02 | expected | 240 |
| A02 | unexpected | 48 |
| A03 | expected | 240 |
| A03 | unexpected | 48 |
kable(head(df.pal %>% count(subID, phase)))
| subID | phase | n |
|---|---|---|
| A01 | postvolatile | 72 |
| A01 | prevolatile | 72 |
| A01 | volatile | 144 |
| A02 | postvolatile | 72 |
| A02 | prevolatile | 72 |
| A02 | volatile | 144 |
kable(head(df.pal %>% count(subID, difficulty)))
| subID | difficulty | n |
|---|---|---|
| A01 | difficult | 96 |
| A01 | easy | 96 |
| A01 | medium | 96 |
| A02 | difficult | 96 |
| A02 | easy | 96 |
| A02 | medium | 96 |
kable(head(df.pal %>% count(subID, expected, phase)))
| subID | expected | phase | n |
|---|---|---|---|
| A01 | expected | postvolatile | 60 |
| A01 | expected | prevolatile | 60 |
| A01 | expected | volatile | 120 |
| A01 | unexpected | postvolatile | 12 |
| A01 | unexpected | prevolatile | 12 |
| A01 | unexpected | volatile | 24 |
kable(head(df.pal %>% count(subID, expected, difficulty)))
| subID | expected | difficulty | n |
|---|---|---|---|
| A01 | expected | difficult | 80 |
| A01 | expected | easy | 80 |
| A01 | expected | medium | 80 |
| A01 | unexpected | difficult | 16 |
| A01 | unexpected | easy | 16 |
| A01 | unexpected | medium | 16 |
kable(head(df.pal %>% count(subID, phase, difficulty)))
| subID | phase | difficulty | n |
|---|---|---|---|
| A01 | postvolatile | difficult | 24 |
| A01 | postvolatile | easy | 24 |
| A01 | postvolatile | medium | 24 |
| A01 | prevolatile | difficult | 24 |
| A01 | prevolatile | easy | 24 |
| A01 | prevolatile | medium | 24 |
kable(head(df.pal %>% count(subID, expected, phase, difficulty)))
| subID | expected | phase | difficulty | n |
|---|---|---|---|---|
| A01 | expected | postvolatile | difficult | 20 |
| A01 | expected | postvolatile | easy | 20 |
| A01 | expected | postvolatile | medium | 20 |
| A01 | expected | prevolatile | difficult | 20 |
| A01 | expected | prevolatile | easy | 20 |
| A01 | expected | prevolatile | medium | 20 |
code = "PAL-rt"
# set the formula
f.pal = brms::bf(rt.cor ~ diagnosis * expected * phase * difficulty +
(expected * phase * difficulty | subID) + (1 | trl))
# informative priors based Lawson et al. and Schad, Betancourt & Vasishth (2019)
priors = c(
prior(normal(6.0, 0.3), class = Intercept),
prior(normal(0.0, 0.5), class = sigma),
prior(normal(0, 0.1), class = sd),
prior(lkj(2), class = cor),
prior(normal(100, 100.0), class = ndt),
prior(normal(0.00, 0.04), class = b)
)
Since this model is so complex and the SBC takes ages to run, we only
run 125 simulations - these already took almost two weeks to run. We ran
them in a separate bash script in smaller batches. Both scripts, R and
bash, are in the helper folder.
if (!(file.exists(file.path(cache_dir, sprintf("df_res_%s.rds", code)))) |
!(file.exists(file.path(cache_dir, sprintf("df_div_%s.rds", code))))) {
# read in the data
dat = readRDS(file = file.path(cache_dir, paste0("dat_", code, ".rds")))
# combine the SBC files to one
ls.res = list.files(path = cache_dir, pattern = sprintf("^res_%s*", code))
for (i in 1:length(ls.res)) {
tmp = readRDS(file.path(cache_dir, ls.res[i]))
if (i == 1) {
res = tmp$result
} else {
res = bind_results(res, tmp$result)
}
}
# put in separate variables
df.results = res$stats
df.backend = res$backend_diagnostics
# save the output
saveRDS(df.results, file.path(cache_dir, sprintf("df_res_%s.rds", code)))
saveRDS(df.backend, file.path(cache_dir, sprintf("df_div_%s.rds", code)))
} else {
# load in the results of the SBC
df.results = readRDS(file.path(cache_dir, sprintf("df_res_%s.rds", code)))
df.backend = readRDS(file.path(cache_dir, sprintf("df_div_%s.rds", code)))
# read in the data
dat = readRDS(file = file.path(cache_dir, paste0("dat_", code, ".rds")))
}
Again, we start by investigating the rhats and the number of divergent samples. This shows that 2 of 125 simulations had at least one parameter that had an rhat of at least 1.05, and 0 models had divergent samples. This suggests that this model performs well and we can continue with our checks by plotting the simulated values to perform prior predictive checks.
# create a matrix out of generated data
dvfakemat = matrix(NA, nrow(dat[['generated']][[1]]), length(dat[['generated']]))
for (i in 1:length(dat[['generated']])) {
dvfakemat[,i] = dat[['generated']][[i]][[1]]
}
# compute one histogram per simulated data-set
dvfakematH = dvfakemat
dvfakematH[dvfakematH > 2000] = 2000
breaks = seq(0, max(dvfakematH, na.rm=T), length.out = 101)
binwidth = breaks[2] - breaks[1]
histmat = matrix(NA, ncol = dim(dvfakemat)[2], nrow = length(breaks)-1)
for (i in 1:dim(dvfakemat)[2]) {
histmat[,i] = hist(dvfakematH[,i], breaks = breaks, plot = F)$counts
}
# for each bin, compute quantiles across histograms
probs = seq(0.1, 0.9, 0.1)
quantmat = as.data.frame(matrix(NA, nrow=dim(histmat)[1], ncol = length(probs)))
names(quantmat) = paste0("p", probs)
for (i in 1:dim(histmat)[1]) {
quantmat[i,] = quantile(histmat[i,], p = probs)
}
quantmat$x = breaks[2:length(breaks)] - binwidth/2 # add bin mean
p1 = ggplot(data = quantmat, aes(x = x)) +
geom_ribbon(aes(ymax = p0.9, ymin = p0.1), fill = c_light) +
geom_ribbon(aes(ymax = p0.8, ymin = p0.2), fill = c_light_highlight) +
geom_ribbon(aes(ymax = p0.7, ymin = p0.3), fill = c_mid) +
geom_ribbon(aes(ymax = p0.6, ymin = p0.4), fill = c_mid_highlight) +
geom_line(aes(y = p0.5), colour = c_dark, linewidth = 1) +
labs(title = "Distribution of simulated values", y = "", x = "") +
theme_bw()
tmpM = apply(dvfakemat, 2, mean) # mean
tmpSD = apply(dvfakemat, 2, sd)
p2 = ggplot() +
stat_bin(aes(x = tmpM), fill = c_dark) +
labs(x = "Mean", title = "Means of simulated values") +
theme_bw()
p3 = ggplot() +
stat_bin(aes(x = tmpSD), fill = c_dark) +
labs(x = "SD", title = "SDs of simulated values") +
theme_bw()
p = ggarrange(p1,
ggarrange(p2, p3, ncol = 2, labels = c("B", "C")),
nrow = 2, labels = "A")
annotate_figure(p, top = text_grob("Prior predictive checks", face = "bold", size = 14))
We are happy with the distribution of the simulated values.
# get simulation numbers with issues
rank = max(df.results$max_rank)
check = merge(df.results %>%
group_by(sim_id) %>%
summarise(
rhat = max(rhat, na.rm = T),
mean_rank = mean(max_rank)
) %>%
filter(rhat >= 1.05 | mean_rank != rank),
df.backend %>% filter(n_divergent > 0), all = T)
# plot SBC with functions from the SBC package focusing on population-level parameters
df.results.b = df.results %>%
filter(substr(variable, 1, 2) == "b_") %>%
filter(!(sim_id %in% check$sim_id))
plot_ecdf_diff(df.results.b) + theme_bw() + theme(legend.position = "none") +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
ggtitle("Computational faithfulness: empirical cumulative distribution function")
plot_rank_hist(df.results.b, bins = 20) + theme_bw() +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
ggtitle("Computational faithfulness: rank histograms")
plot_sim_estimated(df.results.b, alpha = 0.5) + theme_bw() +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
ggtitle("Computational faithfulness: simulated and estimated values")
plot_contraction(df.results.b,
prior_sd = setNames(
c(as.numeric(
gsub(".*, (.+)\\).*", "\\1",
priors[priors$class == "Intercept",]$prior)),
rep(0.04, times = (length(unique(df.results.b$variable))-1))),
unique(df.results.b$variable))) +
theme_bw() +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 5)) +
ggtitle("Computational faithfulness: contraction and z-values")
While some of the ECDFs show discrepancies from the theoretical distribution, we judge all evaluated parameters to be acceptable.
# model formula
f.pup = brms::bf(rel_pupil ~ diagnosis * expected + rts +
(1 | subID)
)
# set weakly informative priors
priors = c(
prior(normal(0, 0.10), class = Intercept),
prior(normal(0.15, 0.15), class = sigma),
prior(normal(0.15, 0.15), class = sd),
prior(normal(0, 0.05), class = b)
)
kable(priors)
| prior | class | coef | group | resp | dpar | nlpar | lb | ub | source |
|---|---|---|---|---|---|---|---|---|---|
| normal(0, 0.1) | Intercept | NA | NA | user | |||||
| normal(0.15, 0.15) | sigma | NA | NA | user | |||||
| normal(0.15, 0.15) | sd | NA | NA | user | |||||
| normal(0, 0.05) | b | NA | NA | user |
code = "PAL-pup-td15"
# extract one dataset from the rt.cor for this
df = dat$generated[[24]] %>%
# add the pupil sizes
mutate(
rel_pupil = rnorm(nrow(dat$generated[[24]]), sd = 0.10),
rts = rt.cor
) %>%
# remove some of the trials to account for data loss
filter(trl <= 252) %>%
select(rel_pupil, diagnosis, expected, subID, rts) %>%
mutate_if(is.character, as.factor)
# print the contrasts
contrasts(df$diagnosis) = contr.sum(2)
contrasts(df$diagnosis)
## [,1]
## A 1
## B -1
contrasts(df$expected) = contr.sum(2)
contrasts(df$expected)
## [,1]
## expected 1
## unexpected -1
# check if the SBC already exists
if (file.exists(file.path(cache_dir, sprintf("df_res_%s.rds", code)))) {
# load in the results of the SBC
df.results = readRDS(file.path(cache_dir, sprintf("df_res_%s.rds", code)))
df.backend = readRDS(file.path(cache_dir, sprintf("df_div_%s.rds", code)))
dat = readRDS(file.path(cache_dir, sprintf("dat_%s.rds", code)))
} else {
# stimulate some data
set.seed(2486)
gen = SBC_generator_brms(f.pup, data = df, prior = priors,
thin = 50, warmup = 10000, refresh = 2000,
generate_lp = TRUE, init = 0.1)
if (file.exists(file.path(cache_dir, sprintf("dat_%s.rds", code)))) {
dat = readRDS(file.path(cache_dir, sprintf("dat_%s.rds", code)))
} else {
dat = generate_datasets(gen, nsim)
saveRDS(dat, file.path(cache_dir, sprintf("dat_%s.rds", code)))
}
# perform the SBC
bck = SBC_backend_brms_from_generator(gen, chains = 4, thin = 1,
warmup = 1500, iter = 6000,
init = 0.1,
control = list(max_treedepth = 15))
res = compute_SBC(dat, bck, keep_fits = F,
cache_mode = "results",
cache_location = file.path(cache_dir, sprintf("res_%s", code)))
df.results = res$stats
df.backend = res$backend_diagnostics
saveRDS(df.results, file = file.path(cache_dir,
paste0("df_res_", code, ".rds")))
saveRDS(df.backend, file = file.path(cache_dir,
paste0("df_div_", code, ".rds")))
}
Again, we start by investigating the rhats and the number of divergent samples. This shows that 134 of 250 simulations had at least one parameter that had an rhat of at least 1.05, and 0 models had divergent samples. This suggests some problems with this model. Let’s investigate this a bit further.
df.results %>% group_by(sim_id) %>% mutate(max_rhat = max(rhat, na.rm = T)) %>% filter(rhat == max_rhat & max_rhat >= 1.05) %>% group_by(variable) %>% count() %>% arrange(desc(n))
## # A tibble: 53 × 2
## # Groups: variable [53]
## variable n
## <chr> <int>
## 1 Intercept 37
## 2 b_diagnosis1 21
## 3 z_1[1,7] 4
## 4 lprior 3
## 5 z_1[1,14] 3
## 6 z_1[1,17] 3
## 7 z_1[1,6] 3
## 8 r_subID[A04,Intercept] 2
## 9 r_subID[A07,Intercept] 2
## 10 r_subID[B01,Intercept] 2
## # ℹ 43 more rows
There are quite some large rhats for the Intercept and the diagnosis effect. We tried increasing the iterations, but it did not help. This might be connected to the maximum treedepth which is saturated in all of this simulations. We will continue for now but watch out for this in the final model.
# create a matrix out of generated data
dvfakemat = matrix(NA, nrow(dat[['generated']][[1]]), length(dat[['generated']]))
for (i in 1:length(dat[['generated']])) {
dvfakemat[,i] = dat[['generated']][[i]][[1]]
}
# compute one histogram per simulated data-set
dvfakematH = dvfakemat
dvfakematH[dvfakematH > 10] = 10
dvfakematH[dvfakematH < -10] = -10
breaks = seq(min(dvfakematH, na.rm=T), max(dvfakematH, na.rm=T), length.out = 101)
binwidth = breaks[2] - breaks[1]
histmat = matrix(NA, ncol = dim(dvfakemat)[2], nrow = length(breaks)-1)
for (i in 1:dim(dvfakemat)[2]) {
histmat[,i] = hist(dvfakematH[,i], breaks = breaks, plot = F)$counts
}
# for each bin, compute quantiles across histograms
probs = seq(0.1, 0.9, 0.1)
quantmat = as.data.frame(matrix(NA, nrow=dim(histmat)[1], ncol = length(probs)))
names(quantmat) = paste0("p", probs)
for (i in 1:dim(histmat)[1]) {
quantmat[i,] = quantile(histmat[i,], p = probs)
}
quantmat$x = breaks[2:length(breaks)] - binwidth/2 # add bin mean
p1 = ggplot(data = quantmat, aes(x = x)) +
geom_ribbon(aes(ymax = p0.9, ymin = p0.1), fill = c_light) +
geom_ribbon(aes(ymax = p0.8, ymin = p0.2), fill = c_light_highlight) +
geom_ribbon(aes(ymax = p0.7, ymin = p0.3), fill = c_mid) +
geom_ribbon(aes(ymax = p0.6, ymin = p0.4), fill = c_mid_highlight) +
geom_line(aes(y = p0.5), colour = c_dark, linewidth = 1) +
labs(title = "Distribution of simulated values", y = "", x = "") +
theme_bw()
tmpM = apply(dvfakemat, 2, mean) # mean
tmpSD = apply(dvfakemat, 2, sd)
p2 = ggplot() +
stat_bin(aes(x = tmpM), fill = c_dark) +
labs(x = "Mean", title = "Means of simulated values") +
theme_bw()
p3 = ggplot() +
stat_bin(aes(x = tmpSD), fill = c_dark) +
labs(x = "SD", title = "SDs of simulated values") +
theme_bw()
p = ggarrange(p1,
ggarrange(p2, p3, ncol = 2, labels = c("B", "C")),
nrow = 2, labels = "A")
annotate_figure(p, top = text_grob("Prior predictive checks", face = "bold", size = 14))
There are quite some extreme values here that are unlikely. However, narrowing the priors also decreased the contraction. Thus, we stick with these rather wide priors.
# get simulation numbers with issues
rank = max(df.results$max_rank)
check = merge(df.results %>%
group_by(sim_id) %>%
summarise(
rhat = max(rhat, na.rm = T),
mean_rank = mean(max_rank)
) %>%
filter(rhat >= 1.05 | mean_rank != rank),
df.backend %>% filter(n_divergent > 0), all = T)
# plot SBC with functions from the SBC package focusing on population-level parameters
df.results.b = df.results %>%
filter(substr(variable, 1, 2) == "b_") %>%
filter(!(sim_id %in% check$sim_id))
plot_ecdf_diff(df.results.b) + theme_bw() + theme(legend.position = "none") +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
ggtitle("Computational faithfulness: empirical cumulative distribution function")
plot_rank_hist(df.results.b, bins = 20) + theme_bw() +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
ggtitle("Computational faithfulness: rank histograms")
plot_sim_estimated(df.results.b, alpha = 0.5) + theme_bw() +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
ggtitle("Computational faithfulness: simulated and estimated values")
plot_contraction(df.results.b,
prior_sd = setNames(
c(as.numeric(
gsub(".*, (.+)\\).*", "\\1",
priors[priors$class == "Intercept",]$prior)),
rep(0.04, times = (length(unique(df.results.b$variable))-1))),
unique(df.results.b$variable))) +
theme_bw() +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 5)) +
ggtitle("Computational faithfulness: contraction and z-values")
Again, we have a few suboptimal values but no bias pattern. Thus, we judge this to be acceptable.
# model formula
f = brms::bf(p ~ diagnosis)
# set weakly informative priors
priors = c(
prior(normal(0, 0.50), class = Intercept),
prior(normal(0, 0.50), class = sigma),
prior(normal(0, 0.25), class = b)
)
kable(priors)
| prior | class | coef | group | resp | dpar | nlpar | lb | ub | source |
|---|---|---|---|---|---|---|---|---|---|
| normal(0, 0.5) | Intercept | NA | NA | user | |||||
| normal(0, 0.5) | sigma | NA | NA | user | |||||
| normal(0, 0.25) | b | NA | NA | user |
# runs faster, so let's do some more simulations
nsim = 1000
code = "PAL-hgf"
# create dataframe
df = df.pal %>%
select(subID, diagnosis) %>%
distinct() %>%
mutate(p = 1)
# print the contrasts
contrasts(df$diagnosis)
## [,1]
## A 1
## B -1
# check if the SBC already exists
if (file.exists(file.path(cache_dir, sprintf("df_res_%s.rds", code)))) {
# load in the results of the SBC
df.results = readRDS(file.path(cache_dir, sprintf("df_res_%s.rds", code)))
df.backend = readRDS(file.path(cache_dir, sprintf("df_div_%s.rds", code)))
dat = readRDS(file.path(cache_dir, sprintf("dat_%s.rds", code)))
} else {
# stimulate some data
set.seed(2486)
gen = SBC_generator_brms(f, data = df, prior = priors,
thin = 50, warmup = 10000, refresh = 2000,
generate_lp = TRUE, init = 0.1)
if (file.exists(file.path(cache_dir, sprintf("dat_%s.rds", code)))) {
dat = readRDS(file.path(cache_dir, sprintf("dat_%s.rds", code)))
} else {
dat = generate_datasets(gen, nsim)
saveRDS(dat, file.path(cache_dir, sprintf("dat_%s.rds", code)))
}
# perform the SBC
bck = SBC_backend_brms_from_generator(gen, chains = 4, thin = 1,
warmup = 1000, iter = 3000)
res = compute_SBC(dat, bck, keep_fits = F,
cache_mode = "results",
cache_location = file.path(cache_dir, sprintf("res_%s", code)))
df.results = res$stats
df.backend = res$backend_diagnostics
saveRDS(df.results, file = file.path(cache_dir,
paste0("df_res_", code, ".rds")))
saveRDS(df.backend, file = file.path(cache_dir,
paste0("df_div_", code, ".rds")))
}
Again, we start by investigating the rhats and the number of divergent samples. This shows that 0 of 1000 simulations had at least one parameter that had an rhat of at least 1.05, and 0 models had divergent samples. This suggests no problems with this model.
# create a matrix out of generated data
dvfakemat = matrix(NA, nrow(dat[['generated']][[1]]), length(dat[['generated']]))
for (i in 1:length(dat[['generated']])) {
dvfakemat[,i] = dat[['generated']][[i]][[1]]
}
# compute one histogram per simulated data-set
dvfakematH = dvfakemat
dvfakematH[dvfakematH > 10] = 10
dvfakematH[dvfakematH < -10] = -10
breaks = seq(min(dvfakematH, na.rm=T), max(dvfakematH, na.rm=T), length.out = 101)
binwidth = breaks[2] - breaks[1]
histmat = matrix(NA, ncol = dim(dvfakemat)[2], nrow = length(breaks)-1)
for (i in 1:dim(dvfakemat)[2]) {
histmat[,i] = hist(dvfakematH[,i], breaks = breaks, plot = F)$counts
}
# for each bin, compute quantiles across histograms
probs = seq(0.1, 0.9, 0.1)
quantmat = as.data.frame(matrix(NA, nrow=dim(histmat)[1], ncol = length(probs)))
names(quantmat) = paste0("p", probs)
for (i in 1:dim(histmat)[1]) {
quantmat[i,] = quantile(histmat[i,], p = probs)
}
quantmat$x = breaks[2:length(breaks)] - binwidth/2 # add bin mean
p1 = ggplot(data = quantmat, aes(x = x)) +
geom_ribbon(aes(ymax = p0.9, ymin = p0.1), fill = c_light) +
geom_ribbon(aes(ymax = p0.8, ymin = p0.2), fill = c_light_highlight) +
geom_ribbon(aes(ymax = p0.7, ymin = p0.3), fill = c_mid) +
geom_ribbon(aes(ymax = p0.6, ymin = p0.4), fill = c_mid_highlight) +
geom_line(aes(y = p0.5), colour = c_dark, linewidth = 1) +
labs(title = "Distribution of simulated values", y = "", x = "") +
theme_bw()
tmpM = apply(dvfakemat, 2, mean) # mean
tmpSD = apply(dvfakemat, 2, sd)
p2 = ggplot() +
stat_bin(aes(x = tmpM), fill = c_dark) +
labs(x = "Mean", title = "Means of simulated values") +
theme_bw()
p3 = ggplot() +
stat_bin(aes(x = tmpSD), fill = c_dark) +
labs(x = "SD", title = "SDs of simulated values") +
theme_bw()
p = ggarrange(p1,
ggarrange(p2, p3, ncol = 2, labels = c("B", "C")),
nrow = 2, labels = "A")
annotate_figure(p, top = text_grob("Prior predictive checks", face = "bold", size = 14))
We are happy with the distribution of the simulated values.
# get simulation numbers with issues
rank = max(df.results$max_rank)
check = merge(df.results %>%
group_by(sim_id) %>%
summarise(
rhat = max(rhat, na.rm = T),
mean_rank = mean(max_rank)
) %>%
filter(rhat >= 1.05 | mean_rank != rank),
df.backend %>% filter(n_divergent > 0), all = T)
# plot SBC with functions from the SBC package focusing on population-level parameters
df.results.b = df.results %>%
filter(substr(variable, 1, 2) == "b_") %>%
filter(!(sim_id %in% check$sim_id))
plot_ecdf_diff(df.results.b) + theme_bw() + theme(legend.position = "none") +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
ggtitle("Computational faithfulness: empirical cumulative distribution function")
plot_rank_hist(df.results.b, bins = 20) + theme_bw() +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
ggtitle("Computational faithfulness: rank histograms")
plot_sim_estimated(df.results.b, alpha = 0.5) + theme_bw() +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
ggtitle("Computational faithfulness: simulated and estimated values")
plot_contraction(df.results.b,
prior_sd = setNames(
c(as.numeric(
gsub(".*, (.+)\\).*", "\\1",
priors[priors$class == "Intercept",]$prior)),
rep(0.25, times = (length(unique(df.results.b$variable))-1))),
unique(df.results.b$variable))) +
theme_bw() +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 5)) +
ggtitle("Computational faithfulness: contraction and z-values")
Again, some deviations but no bias pattern.
# model formula
f = brms::bf(p ~ diagnosis * level * change + (level + change | subID))
# set weakly informative priors
priors = c(
prior(normal(-5, 2.00), class = Intercept),
prior(normal( 0.50, 0.50), class = sigma),
prior(normal( 0, 1.00), class = b),
prior(normal( 0.50, 0.50), class = sd),
prior(lkj(2), class = cor)
)
kable(priors)
| prior | class | coef | group | resp | dpar | nlpar | lb | ub | source |
|---|---|---|---|---|---|---|---|---|---|
| normal(-5, 2) | Intercept | NA | NA | user | |||||
| normal(0.5, 0.5) | sigma | NA | NA | user | |||||
| normal(0, 1) | b | NA | NA | user | |||||
| normal(0.5, 0.5) | sd | NA | NA | user | |||||
| lkj(2) | cor | NA | NA | user |
code = "PAL-lru"
# create dataframe
df = df.pal %>%
select(subID, diagnosis) %>%
distinct() %>%
mutate(p = 1)
df = rbind(df %>% mutate(level = "2", change = "pre2vol"),
df %>% mutate(level = "2", change = "vol2post"),
df %>% mutate(level = "3", change = "pre2vol"),
df %>% mutate(level = "3", change = "vol2post")) %>%
mutate_if(is.character, as.factor)
# print the contrasts
contrasts(df$diagnosis) = contr.sum(2)
contrasts(df$diagnosis)
## [,1]
## A 1
## B -1
contrasts(df$level) = contr.sum(2)
contrasts(df$level)
## [,1]
## 2 1
## 3 -1
contrasts(df$change) = contr.sum(2)
contrasts(df$change)
## [,1]
## pre2vol 1
## vol2post -1
# check if the SBC already exists
if (file.exists(file.path(cache_dir, sprintf("df_res_%s.rds", code)))) {
# load in the results of the SBC
df.results = readRDS(file.path(cache_dir, sprintf("df_res_%s.rds", code)))
df.backend = readRDS(file.path(cache_dir, sprintf("df_div_%s.rds", code)))
dat = readRDS(file.path(cache_dir, sprintf("dat_%s.rds", code)))
} else {
# stimulate some data
set.seed(2486)
gen = SBC_generator_brms(f, data = df, prior = priors, family = lognormal,
thin = 50, warmup = 10000, refresh = 2000,
generate_lp = TRUE, init = 0.1)
if (file.exists(file.path(cache_dir, sprintf("dat_%s.rds", code)))) {
dat = readRDS(file.path(cache_dir, sprintf("dat_%s.rds", code)))
} else {
dat = generate_datasets(gen, nsim)
saveRDS(dat, file.path(cache_dir, sprintf("dat_%s.rds", code)))
}
# perform the SBC
bck = SBC_backend_brms_from_generator(gen, chains = 4, thin = 1,
warmup = 1000, iter = 3000)
res = compute_SBC(dat, bck, keep_fits = F,
cache_mode = "results",
cache_location = file.path(cache_dir, sprintf("res_%s", code)))
df.results = res$stats
df.backend = res$backend_diagnostics
saveRDS(df.results, file = file.path(cache_dir,
paste0("df_res_", code, ".rds")))
saveRDS(df.backend, file = file.path(cache_dir,
paste0("df_div_", code, ".rds")))
}
Again, we start by investigating the rhats and the number of divergent samples. This shows that 5 of 1000 simulations had at least one parameter that had an rhat of at least 1.05, and 21 models had divergent samples. This suggests this model performs fine.
# create a matrix out of generated data
dvfakemat = matrix(NA, nrow(dat[['generated']][[1]]), length(dat[['generated']]))
for (i in 1:length(dat[['generated']])) {
dvfakemat[,i] = dat[['generated']][[i]][[1]]
}
# compute one histogram per simulated data-set
dvfakematH = dvfakemat
dvfakematH[dvfakematH > 1] = 1
#dvfakematH[dvfakematH < -10] = -10
breaks = seq(min(dvfakematH, na.rm=T), max(dvfakematH, na.rm=T), length.out = 101)
binwidth = breaks[2] - breaks[1]
histmat = matrix(NA, ncol = dim(dvfakemat)[2], nrow = length(breaks)-1)
for (i in 1:dim(dvfakemat)[2]) {
histmat[,i] = hist(dvfakematH[,i], breaks = breaks, plot = F)$counts
}
# for each bin, compute quantiles across histograms
probs = seq(0.1, 0.9, 0.1)
quantmat = as.data.frame(matrix(NA, nrow=dim(histmat)[1], ncol = length(probs)))
names(quantmat) = paste0("p", probs)
for (i in 1:dim(histmat)[1]) {
quantmat[i,] = quantile(histmat[i,], p = probs)
}
quantmat$x = breaks[2:length(breaks)] - binwidth/2 # add bin mean
p1 = ggplot(data = quantmat, aes(x = x)) +
geom_ribbon(aes(ymax = p0.9, ymin = p0.1), fill = c_light) +
geom_ribbon(aes(ymax = p0.8, ymin = p0.2), fill = c_light_highlight) +
geom_ribbon(aes(ymax = p0.7, ymin = p0.3), fill = c_mid) +
geom_ribbon(aes(ymax = p0.6, ymin = p0.4), fill = c_mid_highlight) +
geom_line(aes(y = p0.5), colour = c_dark, linewidth = 1) +
labs(title = "Distribution of simulated values", y = "", x = "") +
theme_bw()
tmpM = apply(dvfakemat, 2, mean) # mean
tmpSD = apply(dvfakemat, 2, sd)
p2 = ggplot() +
stat_bin(aes(x = tmpM), fill = c_dark) +
labs(x = "Mean", title = "Means of simulated values") +
theme_bw()
p3 = ggplot() +
stat_bin(aes(x = tmpSD), fill = c_dark) +
labs(x = "SD", title = "SDs of simulated values") +
theme_bw()
p = ggarrange(p1,
ggarrange(p2, p3, ncol = 2, labels = c("B", "C")),
nrow = 2, labels = "A")
annotate_figure(p, top = text_grob("Prior predictive checks", face = "bold", size = 14))
Again, some extreme values but we stick with the wider priors.
# get simulation numbers with issues
rank = max(df.results$max_rank)
check = merge(df.results %>%
group_by(sim_id) %>%
summarise(
rhat = max(rhat, na.rm = T),
mean_rank = mean(max_rank)
) %>%
filter(rhat >= 1.05 | mean_rank != rank),
df.backend %>% filter(n_divergent > 0), all = T)
# plot SBC with functions from the SBC package focusing on population-level parameters
df.results.b = df.results %>%
filter(substr(variable, 1, 2) == "b_") %>%
filter(!(sim_id %in% check$sim_id))
plot_ecdf_diff(df.results.b) + theme_bw() + theme(legend.position = "none") +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
ggtitle("Computational faithfulness: empirical cumulative distribution function")
plot_rank_hist(df.results.b, bins = 20) + theme_bw() +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
ggtitle("Computational faithfulness: rank histograms")
plot_sim_estimated(df.results.b, alpha = 0.5) + theme_bw() +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
ggtitle("Computational faithfulness: simulated and estimated values")
plot_contraction(df.results.b,
prior_sd = setNames(
c(as.numeric(
gsub(".*, (.+)\\).*", "\\1",
priors[priors$class == "Intercept",]$prior)),
rep(1, times = (length(unique(df.results.b$variable))-1))),
unique(df.results.b$variable))) +
theme_bw() +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 5)) +
ggtitle("Computational faithfulness: contraction and z-values")
We judge this to be acceptable.
# model formula
f = brms::bf(p ~ s1 + s1 + s3 + s4 + s5 + s6 + s7 )
# set weakly informative priors
priors = c(
prior(normal(0.5, 0.50), class = Intercept),
prior(normal(0, 1.00), class = b)
)
kable(priors)
| prior | class | coef | group | resp | dpar | nlpar | lb | ub | source |
|---|---|---|---|---|---|---|---|---|---|
| normal(0.5, 0.5) | Intercept | NA | NA | user | |||||
| normal(0, 1) | b | NA | NA | user |
code = "PAL-ber"
# create dataframe
df = df.pal %>%
select(subID, diagnosis) %>%
distinct() %>%
mutate(p = if_else(diagnosis == "A", 1, 0),
s1 = rnorm(length(unique(df.pal$subID))),
s2 = rnorm(length(unique(df.pal$subID))),
s3 = rnorm(length(unique(df.pal$subID))),
s4 = rnorm(length(unique(df.pal$subID))),
s5 = rnorm(length(unique(df.pal$subID))),
s6 = rnorm(length(unique(df.pal$subID))),
s7 = rnorm(length(unique(df.pal$subID)))
)
# check if the SBC already exists
if (file.exists(file.path(cache_dir, sprintf("df_res_%s.rds", code)))) {
# load in the results of the SBC
df.results = readRDS(file.path(cache_dir, sprintf("df_res_%s.rds", code)))
df.backend = readRDS(file.path(cache_dir, sprintf("df_div_%s.rds", code)))
dat = readRDS(file.path(cache_dir, sprintf("dat_%s.rds", code)))
} else {
# stimulate some data
set.seed(2486)
gen = SBC_generator_brms(f, data = df, prior = priors, family = bernoulli,
thin = 50, warmup = 10000, refresh = 2000,
generate_lp = TRUE, init = 0.1)
if (file.exists(file.path(cache_dir, sprintf("dat_%s.rds", code)))) {
dat = readRDS(file.path(cache_dir, sprintf("dat_%s.rds", code)))
} else {
dat = generate_datasets(gen, nsim)
saveRDS(dat, file.path(cache_dir, sprintf("dat_%s.rds", code)))
}
# perform the SBC
bck = SBC_backend_brms_from_generator(gen, chains = 4, thin = 1,
warmup = 1000, iter = 3000)
res = compute_SBC(dat, bck, keep_fits = F,
cache_mode = "results",
cache_location = file.path(cache_dir, sprintf("res_%s", code)))
df.results = res$stats
df.backend = res$backend_diagnostics
saveRDS(df.results, file = file.path(cache_dir,
paste0("df_res_", code, ".rds")))
saveRDS(df.backend, file = file.path(cache_dir,
paste0("df_div_", code, ".rds")))
}
Again, we start by investigating the rhats and the number of divergent samples. This shows that 0 of 250 simulations had at least one parameter that had an rhat of at least 1.05, and 0 models had divergent samples. This suggests this model performs fine.
# create a matrix out of generated data
dvfakemat = matrix(NA, nrow(dat[['generated']][[1]]), length(dat[['generated']]))
for (i in 1:length(dat[['generated']])) {
dvfakemat[,i] = dat[['generated']][[i]][[1]]
}
# compute one histogram per simulated data-set
options = c(0, 1)
histmat = matrix(NA, ncol = dim(dvfakemat)[2], length(options))
for (i in 1:dim(dvfakemat)[2]) {
for (j in 1:length(options))
{
histmat[j,i] = sum(dvfakemat[,i] == options[j])
}
}
# for each bin, compute quantiles across histograms
probs = seq(0.1, 0.9, 0.1)
quantmat= as.data.frame(matrix(NA, nrow=dim(histmat)[1], ncol = length(probs)))
names(quantmat) = paste0("p", probs)
for (i in 1:dim(histmat)[1]) {
quantmat[i,] = quantile(histmat[i,], p = probs)
}
quantmat$x = c("group 0", "group 1")
p1 = ggplot(data = quantmat, aes(x = x)) +
geom_bar(aes(y = p0.9), fill = c_light, stat = "identity") +
geom_bar(aes(y = p0.8), fill = c_light_highlight, stat = "identity") +
geom_bar(aes(y = p0.7), fill = c_mid, stat = "identity") +
geom_bar(aes(y = p0.6), fill = c_mid_highlight, stat = "identity") +
geom_bar(aes(y = p0.5), fill = c_dark, stat = "identity") +
labs(title = "Prior predictive distribution", y = "", x = "") +
theme_bw()
tmpM = apply(dvfakemat, 2, mean) # mean
tmpSD = apply(dvfakemat, 2, sd)
p2 = ggplot() +
stat_bin(aes(x = tmpM), fill = c_dark) +
labs(x = "Mean", title = "Means of simulated values") +
theme_bw()
p3 = ggplot() +
stat_bin(aes(x = tmpSD), fill = c_dark) +
labs(x = "SD", title = "SDs of simulated values") +
theme_bw()
p = ggarrange(p1,
ggarrange(p2, p3, ncol = 2, labels = c("B", "C")),
nrow = 2, labels = "A")
annotate_figure(p, top = text_grob("Prior predictive checks", face = "bold", size = 14))
There is a slight bias for one group, but we judge this to be fine.
# get simulation numbers with issues
rank = max(df.results$max_rank)
check = merge(df.results %>%
group_by(sim_id) %>%
summarise(
rhat = max(rhat, na.rm = T),
mean_rank = mean(max_rank)
) %>%
filter(rhat >= 1.05 | mean_rank != rank),
df.backend %>% filter(n_divergent > 0), all = T)
# plot SBC with functions from the SBC package focusing on population-level parameters
df.results.b = df.results %>%
filter(substr(variable, 1, 2) == "b_") %>%
filter(!(sim_id %in% check$sim_id))
plot_ecdf_diff(df.results.b) + theme_bw() + theme(legend.position = "none") +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
ggtitle("Computational faithfulness: empirical cumulative distribution function")
plot_rank_hist(df.results.b, bins = 20) + theme_bw() +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
ggtitle("Computational faithfulness: rank histograms")
plot_sim_estimated(df.results.b, alpha = 0.5) + theme_bw() +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
ggtitle("Computational faithfulness: simulated and estimated values")
plot_contraction(df.results.b,
prior_sd = setNames(
c(as.numeric(
gsub(".*, (.+)\\).*", "\\1",
priors[priors$class == "Intercept",]$prior)),
rep(1, times = (length(unique(df.results.b$variable))-1))),
unique(df.results.b$variable))) +
theme_bw() +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 5)) +
ggtitle("Computational faithfulness: contraction and z-values")
We judge this to be acceptable.
# model formula
f = brms::bf(p ~ diagnosis * phase + (1 | subID) )
# set weakly informative priors
priors = c(
prior(normal( 2, 1.50), class = Intercept),
prior(normal( 0, 0.50), class = sigma),
prior(normal( 0, 0.50), class = b),
prior(normal( 0, 0.50), class = sd)
)
kable(priors)
| prior | class | coef | group | resp | dpar | nlpar | lb | ub | source |
|---|---|---|---|---|---|---|---|---|---|
| normal(2, 1.5) | Intercept | NA | NA | user | |||||
| normal(0, 0.5) | sigma | NA | NA | user | |||||
| normal(0, 0.5) | b | NA | NA | user | |||||
| normal(0, 0.5) | sd | NA | NA | user |
code = "PAL-ddm"
# create dataframe
df = df.pal %>%
select(subID, diagnosis, phase) %>%
distinct() %>%
mutate(p = 1)
# print the contrasts
contrasts(df$diagnosis) = contr.sum(2)
contrasts(df$diagnosis)
## [,1]
## A 1
## B -1
contrasts(df$phase) = contr.sum(3)
contrasts(df$phase)
## [,1] [,2]
## postvolatile 1 0
## prevolatile 0 1
## volatile -1 -1
# check if the SBC already exists
if (file.exists(file.path(cache_dir, sprintf("df_res_%s.rds", code)))) {
# load in the results of the SBC
df.results = readRDS(file.path(cache_dir, sprintf("df_res_%s.rds", code)))
df.backend = readRDS(file.path(cache_dir, sprintf("df_div_%s.rds", code)))
dat = readRDS(file.path(cache_dir, sprintf("dat_%s.rds", code)))
} else {
# stimulate some data
set.seed(2486)
gen = SBC_generator_brms(f, data = df, prior = priors,
thin = 50, warmup = 10000, refresh = 2000,
generate_lp = TRUE, init = 0.1)
if (file.exists(file.path(cache_dir, sprintf("dat_%s.rds", code)))) {
dat = readRDS(file.path(cache_dir, sprintf("dat_%s.rds", code)))
} else {
dat = generate_datasets(gen, nsim)
saveRDS(dat, file.path(cache_dir, sprintf("dat_%s.rds", code)))
}
# perform the SBC
bck = SBC_backend_brms_from_generator(gen, chains = 4, thin = 1,
warmup = 1000, iter = 3000)
res = compute_SBC(dat, bck, keep_fits = F,
cache_mode = "results",
cache_location = file.path(cache_dir, sprintf("res_%s", code)))
df.results = res$stats
df.backend = res$backend_diagnostics
saveRDS(df.results, file = file.path(cache_dir,
paste0("df_res_", code, ".rds")))
saveRDS(df.backend, file = file.path(cache_dir,
paste0("df_div_", code, ".rds")))
}
Again, we start by investigating the rhats and the number of divergent samples. This shows that 3 of 1000 simulations had at least one parameter that had an rhat of at least 1.05, but 103 models had divergent samples. We will keep this in mind for the final model.
# create a matrix out of generated data
dvfakemat = matrix(NA, nrow(dat[['generated']][[1]]), length(dat[['generated']]))
for (i in 1:length(dat[['generated']])) {
dvfakemat[,i] = dat[['generated']][[i]][[1]]
}
# compute one histogram per simulated data-set
dvfakematH = dvfakemat
#dvfakematH[dvfakematH > 10] = 10
#dvfakematH[dvfakematH < -10] = -10
breaks = seq(min(dvfakematH, na.rm=T), max(dvfakematH, na.rm=T), length.out = 101)
binwidth = breaks[2] - breaks[1]
histmat = matrix(NA, ncol = dim(dvfakemat)[2], nrow = length(breaks)-1)
for (i in 1:dim(dvfakemat)[2]) {
histmat[,i] = hist(dvfakematH[,i], breaks = breaks, plot = F)$counts
}
# for each bin, compute quantiles across histograms
probs = seq(0.1, 0.9, 0.1)
quantmat = as.data.frame(matrix(NA, nrow=dim(histmat)[1], ncol = length(probs)))
names(quantmat) = paste0("p", probs)
for (i in 1:dim(histmat)[1]) {
quantmat[i,] = quantile(histmat[i,], p = probs)
}
quantmat$x = breaks[2:length(breaks)] - binwidth/2 # add bin mean
p1 = ggplot(data = quantmat, aes(x = x)) +
geom_ribbon(aes(ymax = p0.9, ymin = p0.1), fill = c_light) +
geom_ribbon(aes(ymax = p0.8, ymin = p0.2), fill = c_light_highlight) +
geom_ribbon(aes(ymax = p0.7, ymin = p0.3), fill = c_mid) +
geom_ribbon(aes(ymax = p0.6, ymin = p0.4), fill = c_mid_highlight) +
geom_line(aes(y = p0.5), colour = c_dark, linewidth = 1) +
labs(title = "Distribution of simulated values", y = "", x = "") +
theme_bw()
tmpM = apply(dvfakemat, 2, mean) # mean
tmpSD = apply(dvfakemat, 2, sd)
p2 = ggplot() +
stat_bin(aes(x = tmpM), fill = c_dark) +
labs(x = "Mean", title = "Means of simulated values") +
theme_bw()
p3 = ggplot() +
stat_bin(aes(x = tmpSD), fill = c_dark) +
labs(x = "SD", title = "SDs of simulated values") +
theme_bw()
p = ggarrange(p1,
ggarrange(p2, p3, ncol = 2, labels = c("B", "C")),
nrow = 2, labels = "A")
annotate_figure(p, top = text_grob("Prior predictive checks", face = "bold", size = 14))
The distributions look good.
# get simulation numbers with issues
rank = max(df.results$max_rank)
check = merge(df.results %>%
group_by(sim_id) %>%
summarise(
rhat = max(rhat, na.rm = T),
mean_rank = mean(max_rank)
) %>%
filter(rhat >= 1.05 | mean_rank != rank),
df.backend %>% filter(n_divergent > 0), all = T)
# plot SBC with functions from the SBC package focusing on population-level parameters
df.results.b = df.results %>%
filter(substr(variable, 1, 2) == "b_") %>%
filter(!(sim_id %in% check$sim_id))
plot_ecdf_diff(df.results.b) + theme_bw() + theme(legend.position = "none") +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
ggtitle("Computational faithfulness: empirical cumulative distribution function")
plot_rank_hist(df.results.b, bins = 20) + theme_bw() +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
ggtitle("Computational faithfulness: rank histograms")
plot_sim_estimated(df.results.b, alpha = 0.5) + theme_bw() +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 3)) +
ggtitle("Computational faithfulness: simulated and estimated values")
plot_contraction(df.results.b,
prior_sd = setNames(
c(as.numeric(
gsub(".*, (.+)\\).*", "\\1",
priors[priors$class == "Intercept",]$prior)),
rep(0.5, times = (length(unique(df.results.b$variable))-1))),
unique(df.results.b$variable))) +
theme_bw() +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 5)) +
ggtitle("Computational faithfulness: contraction and z-values")
Slight deviations, but no clear bias pattern, thus, fine.